1 Pre-processing

1.1 Load package

library("RColorBrewer")
library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)
#library(SeuratWrappers)
library(harmony)
library(EnsDb.Hsapiens.v86)
library(stringr)
library(dplyr)
library(ggplot2)
library(patchwork)
library(kableExtra)
library(tidyverse)

set.seed(123)

1.2 Parameters

# Paths
path_to_obj <- ("~/Documents/multiome_tonsil_Lucia/results/R_objects/13.tonsil_multiome_bcells_without_doublets_normalized.rds")
path_to_obj_all_data<- ("~/Documents/multiome_tonsil_Lucia/results/R_objects/12.tonsil_multiome_without_cluster7_doublets_normalized.rds")
path_to_markers<-("~/Documents/multiome_tonsil_Lucia/results/tables/tonsil_markers_bcell_075.csv")

1.3 Load data

tonsil_wnn_bcell <- readRDS(path_to_obj)
tonsil_wnn<-readRDS(path_to_obj_all_data)
tonsil_markers_075<-read_csv(path_to_markers)
## New names:
## * `` -> ...1
## Rows: 4688 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): ...1, gene
## dbl (6): p_val, avg_log2FC, pct.1, pct.2, p_val_adj, cluster
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

2 UMAP

 DimPlot(
    tonsil_wnn_bcell,
    group.by = "wsnn_res.0.075",
    reduction = "wnn.umap",
    pt.size = 0.1, label = T
  )

tonsil_wnn_bcell$is_cluster4 <- 
  tonsil_wnn_bcell$wsnn_res.0.075 == "4" 

tonsil_wnn_cluster4 <- subset(tonsil_wnn_bcell, subset = is_cluster4 == TRUE)

barcode_cluster4<-tonsil_wnn_cluster4@meta.data$lib_name_barcode

3 B cell UMAP

3.1 WNN

DimPlot(object = tonsil_wnn_bcell, cells.highlight = barcode_cluster4, reduction = "wnn.umap", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()

3.2 ATAC

DimPlot(object = tonsil_wnn_bcell, cells.highlight = barcode_cluster4, reduction = "umap.atac", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()

3.3 RNA

DimPlot(object = tonsil_wnn_bcell, cells.highlight = barcode_cluster4, reduction = "umap.rna", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()

3.4 split

joint.umap_age<- DimPlot(tonsil_wnn_bcell, label = FALSE, split.by = "age_group", pt.size = 0.1,  reduction = "umap.rna") + plot_annotation(title = 'RNA UMAP ')+ ggtitle('age group') 

joint.umap_hospital<- DimPlot(tonsil_wnn_bcell, label = FALSE, split.by = "hospital", pt.size = 0.1,  reduction = "umap.rna") + plot_annotation(title = 'RNA UMAP ')+ ggtitle('hospital ') 

joint.umap_age

joint.umap_hospital

joint.umap_age<- DimPlot(tonsil_wnn_bcell, label = FALSE, split.by = "age_group", pt.size = 0.1,  reduction = "umap.atac") + plot_annotation(title = 'ATAC UMAP ')+ ggtitle('Age group') 

joint.umap_hospital<- DimPlot(tonsil_wnn_bcell, label = FALSE, split.by = "hospital", pt.size = 0.1,  reduction = "umap.atac") + plot_annotation(title = 'ATAC UMAP ')+ ggtitle('Hospital ') 

joint.umap_age

joint.umap_hospital

4 All data UMAP (without doublets)

4.1 WNN

DimPlot(object = tonsil_wnn, cells.highlight = barcode_cluster4, reduction = "wnn.umap", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()

4.2 ATAC

DimPlot(object = tonsil_wnn, cells.highlight = barcode_cluster4, reduction = "umap.atac", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()

4.3 RNA

DimPlot(object = tonsil_wnn, cells.highlight = barcode_cluster4, reduction = "umap.rna", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()

4.4 split

joint.umap_age<- DimPlot(tonsil_wnn, label = FALSE, split.by = "age_group", pt.size = 0.1,  reduction = "umap.rna") + plot_annotation(title = 'RNA UMAP - All data ')+ ggtitle('age group') 

joint.umap_hospital<- DimPlot(tonsil_wnn, label = FALSE, split.by = "hospital", pt.size = 0.1,  reduction = "umap.rna") + plot_annotation(title = 'RNA UMAP ')+ ggtitle('hospital ') 

joint.umap_age

joint.umap_hospital

joint.umap_age<- DimPlot(tonsil_wnn, label = FALSE, split.by = "age_group", pt.size = 0.1,  reduction = "umap.atac") + plot_annotation(title = 'ATAC UMAP ')+ ggtitle('Age group') 

joint.umap_hospital<- DimPlot(tonsil_wnn, label = FALSE, split.by = "hospital", pt.size = 0.1,  reduction = "umap.atac") + plot_annotation(title = 'ATAC UMAP - All data ')+ ggtitle('Hospital ') 

joint.umap_age

joint.umap_hospital

5 Markers

top10_tonsil_markers_075<-tonsil_markers_075 %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

top10mark_cluster4<-top10_tonsil_markers_075[["gene"]][41:50]
top7_tonsil_markers_075<-tonsil_markers_075 %>% group_by(cluster) %>% top_n(n = 7, wt = avg_log2FC)

top7mark_cluster4<-top7_tonsil_markers_075[["gene"]][29:35]
df_top10<-as.data.frame(top10_tonsil_markers_075)

DT::datatable(df_top10)
#install.packages("htmlwidgets", type = "binary")
#install.packages("DT", type = "binary")

df_top10<-as.data.frame(top10mark_cluster4)

df_top10
##    top10mark_cluster4
## 1                FYB1
## 2              INPP4B
## 3              THEMIS
## 4               PRKCH
## 5                IL7R
## 6                LEF1
## 7                TC2N
## 8               CAMK4
## 9              BCL11B
## 10                ITK
#DoHeatmap(tonsil_wnn_bcell)


DoHeatmap(tonsil_wnn_bcell,top7_tonsil_markers_075$gene,assay = "RNA")
## Warning in DoHeatmap(tonsil_wnn_bcell, top7_tonsil_markers_075$gene, assay
## = "RNA"): The following features were omitted as they were not found in the
## scale.data slot for the RNA assay: FCER2, ZDHHC14

markers_kn<-c("BANK1", "FCER2","CD3D", "IL7R","MKI67", "TOP2A","MARCKSL1", "RGS13", "LMO2", "CCDC88A","GNLY", "NKG7","CD8A",    "FCRL4", "FCRL5", "PLAC8", "SOX5","IGHG1", "IGHA1", "XBP1","LYZ", "S100A8", "CD3D")
DoHeatmap(tonsil_wnn_bcell,markers_kn,assay = "RNA")
## Warning in DoHeatmap(tonsil_wnn_bcell, markers_kn, assay = "RNA"): The following
## features were omitted as they were not found in the scale.data slot for the RNA
## assay: S100A8, PLAC8, CD8A, MARCKSL1, FCER2, BANK1

We are going to represent the top 7 markers

markers_gg_bcell <- function(x,umap){purrr::map(x, function(x) {
  
  p <- FeaturePlot(
    tonsil_wnn_bcell,
    features = x,
    reduction = umap,
    pt.size = 0.1
  )
  p
})}
markers_gg_alldata <- function(x,umap){purrr::map(x, function(x) {
  
  p <- FeaturePlot(
    tonsil_wnn,
    features = x,
    reduction = umap,
    pt.size = 0.1
  )
  p
})}

5.1 Pauli markers

naive_mem_bcell<-c("BANK1", "FCER2")
cd4_tcell<-c("CD3D", "IL7R")
dz_gc_bcell<-c("MKI67", "TOP2A")
lz_gc_bcell<-c("MARCKSL1", "RGS13", "LMO2", "CCDC88A")  
cytotoxic<-c("GNLY", "NKG7", "GZMK", "CD8A")
memory_bcell<-c(    "FCRL4", "FCRL5", "PLAC8", "SOX5")
pc<-c("IGHG1", "IGHA1", "JCHAIN", "XBP1")
myeloid<-c("LYZ", "S100A8") 
poor_q_doublets <-c("FDCSP", "CLU", "CXCL13", "CR2")
doublet_proliferative_tcell<-c("MT2A", "CD3D", "TRAC", "PCNA")
Unk<-c("PTPRCAP", "CD37", "CD74")   
PDC<-c("PTCRA", "LILRA4", "IRF7")

5.2 All data

5.2.1 WNN UMAP

DimPlot(
  tonsil_wnn,
  reduction = "wnn.umap",
  pt.size = 0.1
) + ggtitle('WNN UMAP')

markers_gg_alldata(top7mark_cluster4,"wnn.umap")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

5.2.2 ATAC UMAP

DimPlot(
  tonsil_wnn,
  reduction = "umap.atac",
  pt.size = 0.1
) + ggtitle('ATAC UMAP Harmony')

markers_gg_alldata(top7mark_cluster4,"umap.atac")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

markers_gg_alldata(c("TCL1A","IGHD","FCER2","IGHM","IL4R"),"umap.atac")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

markers_gg_alldata(naive_mem_bcell,"umap.atac")
## [[1]]

## 
## [[2]]

markers_gg_alldata(memory_bcell,"umap.atac")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

markers_gg_alldata(dz_gc_bcell,"umap.atac")
## [[1]]

## 
## [[2]]

markers_gg_alldata(lz_gc_bcell,"umap.atac")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

markers_gg_alldata(pc,"umap.atac")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

5.2.2.1 T cell

markers_gg_alldata(cytotoxic,"umap.atac")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

markers_gg_alldata(cd4_tcell,"umap.atac")
## [[1]]

## 
## [[2]]

markers_gg_alldata(myeloid,"umap.atac")
## [[1]]

## 
## [[2]]

5.2.3 RNA UMAP

DimPlot(
  tonsil_wnn,
  reduction = "umap.rna",
  pt.size = 0.1
) + ggtitle('RNA Harmony')

markers_gg_alldata(top7mark_cluster4,"umap.rna")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

markers_gg_alldata(c("TCL1A","IGHD","FCER2","IGHM","IL4R"),"umap.rna")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

markers_gg_alldata(naive_mem_bcell,"umap.rna")
## [[1]]

## 
## [[2]]

markers_gg_alldata(memory_bcell,"umap.rna")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

markers_gg_alldata(dz_gc_bcell,"umap.rna")
## [[1]]

## 
## [[2]]

markers_gg_alldata(lz_gc_bcell,"umap.rna")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

markers_gg_alldata(pc,"umap.rna")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

5.2.3.1 T cell

markers_gg_alldata(cytotoxic,"umap.rna")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

markers_gg_alldata(cd4_tcell,"umap.rna")
## [[1]]

## 
## [[2]]

markers_gg_alldata(myeloid,"umap.rna")
## [[1]]

## 
## [[2]]

5.3 B cell data

5.3.1 WNN UMAP

markers_gg_bcell(top7mark_cluster4,"wnn.umap")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

5.3.2 ATAC UMAP

DimPlot(
  tonsil_wnn_bcell,
  reduction = "umap.atac",
  pt.size = 0.1
) + ggtitle('ATAC Harmony B cell')

markers_gg_bcell(top7mark_cluster4,"umap.atac")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

5.3.3 RNA UMAP

DimPlot(
  tonsil_wnn_bcell,
  reduction = "umap.rna",
  pt.size = 0.1
) + ggtitle('RNA Harmony B cell')

markers_gg_bcell(top7mark_cluster4,"umap.rna")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

6 Session Information

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] es_ES.UTF-8/es_ES.UTF-8/es_ES.UTF-8/C/es_ES.UTF-8/es_ES.UTF-8
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] forcats_0.5.1             purrr_0.3.4              
##  [3] readr_2.1.2               tidyr_1.2.0              
##  [5] tibble_3.1.6              tidyverse_1.3.1          
##  [7] kableExtra_1.3.4          patchwork_1.1.1          
##  [9] ggplot2_3.3.5             dplyr_1.0.8              
## [11] stringr_1.4.0             EnsDb.Hsapiens.v86_2.99.0
## [13] ensembldb_2.18.2          AnnotationFilter_1.18.0  
## [15] GenomicFeatures_1.46.1    AnnotationDbi_1.56.2     
## [17] Biobase_2.54.0            harmony_0.1.0            
## [19] Rcpp_1.0.8                future_1.23.0            
## [21] GenomicRanges_1.46.1      GenomeInfoDb_1.30.0      
## [23] IRanges_2.28.0            S4Vectors_0.32.3         
## [25] BiocGenerics_0.40.0       SeuratObject_4.0.4       
## [27] Seurat_4.1.0              Signac_1.5.0             
## [29] RColorBrewer_1.1-2        BiocStyle_2.22.0         
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3              SnowballC_0.7.0            
##   [3] rtracklayer_1.54.0          scattermore_0.7            
##   [5] bit64_4.0.5                 knitr_1.37                 
##   [7] irlba_2.3.5                 DelayedArray_0.20.0        
##   [9] data.table_1.14.2           rpart_4.1-15               
##  [11] KEGGREST_1.34.0             RCurl_1.98-1.5             
##  [13] generics_0.1.2              cowplot_1.1.1              
##  [15] RSQLite_2.2.9               RANN_2.6.1                 
##  [17] bit_4.0.4                   tzdb_0.2.0                 
##  [19] spatstat.data_2.1-2         webshot_0.5.2              
##  [21] xml2_1.3.3                  lubridate_1.8.0            
##  [23] httpuv_1.6.5                SummarizedExperiment_1.24.0
##  [25] assertthat_0.2.1            xfun_0.29                  
##  [27] hms_1.1.1                   jquerylib_0.1.4            
##  [29] evaluate_0.14               promises_1.2.0.1           
##  [31] fansi_1.0.2                 restfulr_0.0.13            
##  [33] progress_1.2.2              dbplyr_2.1.1               
##  [35] readxl_1.3.1                igraph_1.2.11              
##  [37] DBI_1.1.2                   htmlwidgets_1.5.4          
##  [39] sparsesvd_0.2               spatstat.geom_2.3-1        
##  [41] ellipsis_0.3.2              crosstalk_1.2.0            
##  [43] backports_1.4.1             bookdown_0.24              
##  [45] biomaRt_2.50.1              deldir_1.0-6               
##  [47] MatrixGenerics_1.6.0        vctrs_0.3.8                
##  [49] ROCR_1.0-11                 abind_1.4-5                
##  [51] cachem_1.0.6                withr_2.4.3                
##  [53] ggforce_0.3.3               vroom_1.5.7                
##  [55] sctransform_0.3.3           GenomicAlignments_1.30.0   
##  [57] prettyunits_1.1.1           goftest_1.2-3              
##  [59] svglite_2.0.0               cluster_2.1.2              
##  [61] lazyeval_0.2.2              crayon_1.5.0               
##  [63] pkgconfig_2.0.3             slam_0.1-49                
##  [65] labeling_0.4.2              tweenr_1.0.2               
##  [67] nlme_3.1-153                ProtGenerics_1.26.0        
##  [69] rlang_1.0.1                 globals_0.14.0             
##  [71] lifecycle_1.0.1             miniUI_0.1.1.1             
##  [73] filelock_1.0.2              BiocFileCache_2.2.0        
##  [75] modelr_0.1.8                cellranger_1.1.0           
##  [77] polyclip_1.10-0             matrixStats_0.61.0         
##  [79] lmtest_0.9-39               Matrix_1.3-4               
##  [81] ggseqlogo_0.1               zoo_1.8-9                  
##  [83] reprex_2.0.1                ggridges_0.5.3             
##  [85] png_0.1-7                   viridisLite_0.4.0          
##  [87] rjson_0.2.20                bitops_1.0-7               
##  [89] KernSmooth_2.23-20          Biostrings_2.62.0          
##  [91] blob_1.2.2                  parallelly_1.30.0          
##  [93] scales_1.1.1                memoise_2.0.1              
##  [95] magrittr_2.0.2              plyr_1.8.6                 
##  [97] ica_1.0-2                   zlibbioc_1.40.0            
##  [99] compiler_4.1.2              BiocIO_1.4.0               
## [101] fitdistrplus_1.1-6          Rsamtools_2.10.0           
## [103] cli_3.1.1                   XVector_0.34.0             
## [105] listenv_0.8.0               pbapply_1.5-0              
## [107] MASS_7.3-54                 mgcv_1.8-38                
## [109] tidyselect_1.1.1            stringi_1.7.6              
## [111] highr_0.9                   yaml_2.2.2                 
## [113] ggrepel_0.9.1               grid_4.1.2                 
## [115] sass_0.4.0                  fastmatch_1.1-3            
## [117] tools_4.1.2                 future.apply_1.8.1         
## [119] parallel_4.1.2              rstudioapi_0.13            
## [121] lsa_0.73.2                  gridExtra_2.3              
## [123] farver_2.1.0                Rtsne_0.15                 
## [125] digest_0.6.29               BiocManager_1.30.16        
## [127] shiny_1.7.1                 qlcMatrix_0.9.7            
## [129] broom_0.7.12                later_1.3.0                
## [131] RcppAnnoy_0.0.19            httr_1.4.2                 
## [133] colorspace_2.0-2            rvest_1.0.2                
## [135] XML_3.99-0.8                fs_1.5.2                   
## [137] tensor_1.5                  reticulate_1.24            
## [139] splines_4.1.2               uwot_0.1.11                
## [141] RcppRoll_0.3.0              spatstat.utils_2.3-0       
## [143] plotly_4.10.0               systemfonts_1.0.3          
## [145] xtable_1.8-4                jsonlite_1.7.3             
## [147] R6_2.5.1                    pillar_1.7.0               
## [149] htmltools_0.5.2             mime_0.12                  
## [151] glue_1.6.1                  fastmap_1.1.0              
## [153] DT_0.20                     BiocParallel_1.28.3        
## [155] codetools_0.2-18            utf8_1.2.2                 
## [157] lattice_0.20-45             bslib_0.3.1                
## [159] spatstat.sparse_2.1-0       curl_4.3.2                 
## [161] leiden_0.3.9                magick_2.7.3               
## [163] survival_3.2-13             rmarkdown_2.11             
## [165] docopt_0.7.1                munsell_0.5.0              
## [167] GenomeInfoDbData_1.2.7      haven_2.4.3                
## [169] reshape2_1.4.4              gtable_0.3.0               
## [171] spatstat.core_2.3-2